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ABSTRACT 

Feedback from massive stars is believed to play a critical role in driving galactic super- 
winds that enrich the intergalactic medium and shape the galaxy mass function, mass- 
metallicity relation, and other global galaxy properties. In previous papers, we have intro- 
duced new numerical methods for implementing stellar feedback on sub-GMC through galac- 
tic scales in numerical simulations of galaxies; the key physical processes include radiation 
pressure in the UV through IR, supernovae (Type-I & II), stellar winds ("fast" O star through 
"slow" AGB winds), and HII photoionization. Here, we show that these feedback mechanisms 
drive galactic winds with outflow rates as high as ~ 10 — 20 times the galaxy star formation 
rate. The mass-loading efficiency (wind mass loss rate divided by the star formation rate) 
scales roughly as M^ind/^* ^ (where Vc is the galaxy circular velocity), consistent with 
simple momentum-conservation expectations. We use our suite of simulations to study the rel- 
ative contribution of each feedback mechanism to the generation of galactic winds in a range 
of galaxy models, from SMC-like dwarfs and Milky- way analogues to z ~ 2 clumpy disks. In 
massive, gas-rich systems (local starbursts and high-z galaxies), radiation pressure dominates 
the wind generation. By contrast, for MW-like spirals and dwarf galaxies the gas densities 
are much lower and sources of shock-heated gas such as supernovae and stellar winds dom- 
inate the production of large-scale outflows. In all of our models, however, the winds have 
a complex multi-phase structure that depends on the interaction between multiple feedback 
mechanisms operating on different spatial and time scales: any single feedback mechanism 
fails to reproduce the winds observed. We use our simulations to provide fitting functions to 
the wind mass-loading and velocities as a function of galaxy properties, for use in cosmo- 
logical simulations and semi-analytic models. These differ from typically-adopted formulae 
with an explicit dependence on the gas surface density that can be very important in both 
low-density dwarf galaxies and high-density gas -rich galaxies. 

Key words: galaxies: formation — star formation: general — galaxies: evolution — galaxies: 
active — cosmology: theory 



1 INTRODUCTION 

Feedback from massive stars is critical to the evolution of galaxies. 
In cosmological models of galaxy evolution without strong stellar 
feedback, gas rapidly cools and turns into stars, leading to galax- 
ies with star formation rates much higher than observed, and ~ ten 
times the stellar mass found in real galaxies (e.g. !Katz et a l.|1996[ 
Some rville & Primack 19991 [Cole eraLl[2050l fSpringel & Hern-| 
quist||2003b[ |Keres et a l. 2009a and references therein). Simply 
suppressing the rate of star formation does not solve the problem: 
the amount of baryons in real galactic disks is much lower than the 
amount of cool gas in disks found in cosmological simulations, es- 
pecially in low-mass galaxies ( [White & Frenk|199T] for a recent re- 
view see Keres et al. 2009b ). Constraints from the mass-metallicity 
relation and enrichment of the IGM also imply that the baryons 
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cannot simply be prevented from entering galaxy halos along with 
dark matter (Tremonti et al. 2004, Erb et al.||2006[ [Aguirre et aT] 
2001, Pettini et al. 2003, Songaila 2005J . Some process must very 
efficiently remove baryons from galaxies. 

Related problems appear on smaller spatial and time scales. 
The Kennicutt-Schmidt (KS) law implies that star formation is very 
slow within galaxies, with a gas consumption time of ^ 50 dynam- 
ical times ( Ken nicutt||1998 ). Moreover, the integrated fraction of 
mass turned into stars in GMCs over their lifetime is only a few 
to several percent (Zuckerman & Evans 1974, Willi ams & McKee| 
p^97l|Evans^ 1999,, Evans e t al. 2009). Without strong stellar feed- 
back, however, self-gravitating collapse leads to most of the gas 
turning into stars in just a few dynamical times. 

The problem, then, on both galactic and sub-galactic scales, 
is twofold. First, star formation must be "slowed down" at a given 
global/local gas surface density. But this alone would still violate 
integral constraints, producing galaxies and star clusters more mas- 
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sive than observed by an order of magnitude. Thus the second prob- 
lem: on small scales gas must be expelled from GMCs, and on 
galactic scales either prevented from entering, or, more likely in 
our judgment, removed from, the host galaxy. In other words, local 
outflows and global super- winds must be generated that can remove 
gas at a rate rapid compared to the star formation rate. 

Because low-mass galaxies are preferentially baryon-poor, 
matching the faint end of the observed galaxy mass function in cos- 
mological simulations requires that the global efficiency of galactic 
super-winds scales as a declining power of galaxy mass or circu- 
lar velocity. O ppenheimer & Davll ( |2006| ) and |Oppenheimer et al.| 
([2010 ) find that an average scaling Mwind/M* oc V~ produces rea- 
sonably good agreement with the observed mass functions at differ- 
ent redshiftsQwith a normalization such that an SMC-mass dwarf 
has a mass loading Mwind/M* ^10 (although this scaling may still 
over-produce the number of very low-mass galaxies). Large mass- 
loading factors of several times the SFR are also estimated in even 
relatively massive local galaxies and massive star-forming regions 
at z - 2 - 3 ( p^t in 1999 2006, Heckma n et aL|2000 , Newman 
et al.|2012l|Sato et al.|2009nChen et al.|20T0l[Steideret al.|2010| 
Coil et al.|20n] ). 

The scaling Mwind/M* oc is expected from momentum 
conservation arguments, given a number of simplifying assump- 
tions and sufficient global momentum input from supernovae, stel- 
lar winds, radiation pressure, etc. ( [Murray et al.pOOS] ). Direct ob- 
servations, while uncertain, tend to favor velocity and Mwind scal- 
ings similar to this constraint for the bulk of the outflowing gas 
( |Martin|2005]|Rupke et al.|2005||Weiner et al.|2009| ). 

To date, however, numerical simulations have generally not 
been able to produce, from an a priori model, winds with either such 
large absolute mass loading factors or the scaling of mass-loading 
with galaxy mass/velocity. Many simulations, lacking the ability to 
directly resolve the relevant feedback processes, put in winds "by 
hand" by e.g. forcing an outflow rate that scales in a user-specified 
manner with the star formation rate or other par ameters (|Springel| 
& Hernquist|2003a[[Oppenheimer & Dave|2008[ [Sales et aLpoTof 
Genel et al.|2012| ). Alternatively, models that self-consistently in- 
clude stellar feedback have generally been limited to a small subset 
of the relevant processes; the vast majority include only thermal 
feedback via supernovae (i.e. thermal energy injection with some 
average rate that scales with the mass in young stars). However, 
thermal feedback is very inefficient in the dense regions where star 
formation occurs, and in the ISM more broadly in gas-rich galaxies. 
For this reason, such models require further changes to the physics 
in order for thermal energy injection to have a significant effect. 
Often cooling (along with star formation and other hydrodynamic 
processes) is "turned off" for an extended period of time ([Thacker| 
|& Couchman|2000[[Govemato et al.|2007| [Brook et al. ^201 1 ). With 
or without these adjustments, however, such models generally ob- 
tain winds that are weaker than those required to explain the galaxy 
mass function, especially at low mass es (see e.g. [ Guo et al.|2010[ 
[Powell et al.pOTT] [Brook et al.|2011|[Nagamine|2010| and refer- 
ences therein). 

In our view, part of the resolution of this difficulty lies in the 
treatment of the ISM physics within galaxies. Feedback processes 
other than supernovae are critical for suppressing star formation in 



^ Note that the agreement we speak of is with the faint sub-L* end of the 
mass function. It is widely agreed that different physics, perhaps AGN feed- 
back, is critical for the regulation of the bright super-L* end of the mass 
function. We focus here on stellar feedback and star-forming systems. 



dense gas; these include protostellar jets, HII regions, stellar winds, 
and radiation pressure from young stars. Including these mecha- 
nisms self-consistently maintains a reasonable fraction of the ISM 
at densities where the thermal heating from supernovae has a larger 
effect. This conclusion implies that (not surprisingly) a realistic 
treatment of galactic winds requires a more realistic treatment of 
the stellar feedback processes that maintain the multi-phase struc- 
ture of the ISM of galaxies. 

Motivated by this perspective, in [Hopkins et al.[ ( [201 1| ) (Pa- 
per I) and Hopk ins et al.[p012] ) (Paper II) we developed a new set 
of numerical models to follow feedback on small scales in GMCs 
and star-forming regions, in simulations with pc-scale resolution]^ 
These simulations include the momentum imparted locally (on sub- 
GMC scales) from stellar radiation pressure, radiation pressure on 
larger scales via the light that escapes star-forming regions, HII 
photoionization heating, as well as the heating, momentum deposi- 
tion, and mass loss by SNe (Type-I and Type-II) and stellar winds 
(O star and AGB). The feedback is tied to the young stars, with the 
energetics and time-dependence taken directly from stellar evolu- 
tion models. Our models also include realistic cooling to tempera- 
tures < lOOK, and a treatment of the molecular/atomic transition in 
gas and its effect on star formation (following [Krumholz & Gnedin[ 
[20TT] ). 

We showed in Papers I & II that these feedback mechanisms 
produce a quasi- steady ISM in which giant molecular clouds form 
and disperse rapidly, after turning just a few percent of their mass 
into stars. This leads to an ISM with phase structure, turbulent ve- 
locity dispersions, scale heights, and GMC properties (mass func- 
tions, sizes, scaling laws) in reasonable agreement with observa- 
tions. In this paper, we use these same models of stellar feedback to 
quantitatively predict the elusive winds invoked in almost all galaxy 
formation models. 

The remainder of this paper is organized as follows. In ^we 
summarize the galaxy models we use and our methods of imple- 
menting stellar feedback. In § [3] we discuss how each feedback 
mechanism affects the morphology, phase structure, and velocity 
distribution of galactic winds. In §|4] we discuss the mass-loading 
of winds and how the mass loss rate depends on the inclusion of dif- 
ferent feedback mechanisms. We further determine how the outflow 
rate scales with galaxy properties and use our simulations to derive 
more accurate approximations to wind scalings for use in cosmo- 
logical simulations and semi-analytic models. In §[5]we summarize 
our results and discuss their implications. 

2 METHODS 

The simulations used here are described in detail in Paper I (see 
their Section 2 and Tables 1-3) and Paper II (their Section 2). How- 
ever we briefly summarize the most important properties of the 
models here. The simulations were performed with the parallel 
TreeSPH code GADGET-3 ( Springel 2005 ). They include stars, dark 
matter, and gas, with cooling, star formation, and stellar feedback. 

Gas follows a standard atomic cooling curve but in addition 
can cool to < 100 K via fine-structure cooling. This allows it to 
collapse to very high densities, and star formation occurs in dense 
regions above a threshold n > lOOOcm"^, with a rate = ep/fe 
where is the free-fall time and e = 1.5% is an efficiency taken 
from observations of star-forming regions with the same densities 
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( [Krumholz & Tan|2007| and references therein). We further follow 
[Krumholz & Gnedin] ( |2011| ) and calculate the molecular fraction 
within the dense gas as a function of the local column density and 
metallicity, and allow star formation only from that gas. In Paper 
I we show that the SFR is essentially independent of the small- 
scale star formation law (robust to large variations in the thresh- 
old and efficiency, and changes in the power-law index p oc p^~^). 
This is because star formation is feedback-regulated and dense star- 
forming regions grow in mass until sufficient new stars have formed 
to halt further collapse. Likewise in Paper II we show that the cor- 
rections from the molecular chemistry are negligible at the masses 
and metallicities we model. 

2.1 Disk Models 

Our calculations span four distinct initial disk models, designed to 
represent a range of characteristic galaxy types. Each initial disk 
has a bulge, stellar and gaseous disk, and dark matter halo. The 
disks are initialized in equilibrium so that in the absence of cool- 
ing, star formation, and feedback there are no significant transients. 
The gaseous disk is initially vertically pressure- supported, but this 
thermal energy is radiated away in much less than a dynamical time 
and the emergent vertical structure depends on feedback. Our "low" 
resolution runs (used to evolve the simulations for several Gyr, to 
ensure steady-state behavior) use ^ 3 x 10^ particles, with ^ 
particles in the disk, giving SPH smoothing lengths of ^ lOpc in 
the central few kpc of a MW-like disk (the smoothing length scales 
linearly with the disk size/mass scale). Our "standard" resolution 
cases use ~ 30 times as many particles, and correspondingly have 
~ 1 — 5pc smoothing lengths and particle masses of 500 Mq; these 
are run for a few orbital times each. A few ultra-high resolution runs 
used for convergence tests employ ^10^ particles, with sub-pc res- 
olution on kpc scales. 

(1) SMC: an SMC-like dwarf, with baryonic mass Mbar = 
8.9 X lO^M0_andhalomass Mhaio = 2 x 10^^ Mq (concentration 
c= 15), a |Hernquist ( 1990| ) profile bulge with a mass m/, — 10^ Mq, 
and exponential stellar (rrid = 1.3 x IO^Mq) and gas disks (m^ = 
7.5 X IO^Mq) with scale-lengths hd = 0.1 and hg = 2.1 kpc, re- 
spectively. The initial stellar scale-height is zo = 140 pc and the stel- 
lar disk is initialized such that the Toomre 2=1 everywhere. The 
gas and stars are initialized with uniform metallicity Z = 0.1 Zq. 

(2) MW: a MW-like galaxy, with halo and bary- 
onic properties of (Mhaio, c) = (1.6 x 10^^ Mq, 12) and 
(Mbar, mi,, m J, m J = (7.1,1.5,4.7,0.9) x 10^^ Mq, Z = Zq, 
and scale-lengths {hd, hg, zo) = (3.0, 6.0, 0.3) kpc. 

(3) Sbc: a LIRG-like galaxy (i.e. a more gas-rich spiral than is 
characteristic of those observed at low redshifts) with (Mhaio, c) = 
(1.5 X IO^^Mq, 11), {Mb^r,mb,md,mg) (10.5, 1.0,4.0,5.5) x 
lO'^ Mq, Z = 0.3Zq, and {hd, hg, zo) = (1.3, 2.6, 0.13) kpc. 

(4) HiZ: a high-redshift massive starburst disk, chosen to 
match the properties of the observed non-merging but rapidly star- 
forming SMG population, with (Mhaio, c) = (1.4 x IO^^Mq, 3.5) 
and a virial radius appropriately rescaled for a halo at z = 2 rather 
than z = 0, (Mbar, mt, md, nig) (10.7, 0.7, 3, 7) x IO^^Mq, Z 
0.5Zq, and {hd, hg, zo) = (1.6, 3.2, 0.32) kpc. 

2.2 Feedback Models 

The most important physics in these simulations is the model of 
stellar feedback. We include feedback from a variety of mecha- 
nisms, each of which we briefly describe below. More details about 
our implementations of this physics are given in Paper I and Paper 
II. We use a Kroupa (2002 ) initial mass function (IMF) through- 
out and use STARBURST99 ( Leitherer etal.,1999j to calculate the 



stellar luminosity, mass return from stellar winds, supernova rate, 
etc. as a function of the age and metallicity of each star particle. 

(1) Local Momentum Deposition from Radiation Pressure, 
Supernovae, & Stellar Winds: In Paper I, we present the radiation 
pressure aspect of this model for feedback from young star clusters 
in detail. At each timestep, gas particles identify the nearest density 
peak representing the center of the nearest star-forming "clump" or 
GMC-analog. We calculate the total luminosity of the star particles 
inside the sphere defined by the distance from the center of this 
star-forming region to the gas particle of interest; the incident flux 
on the gas is then determined assuming that the local star forming 
region is optically thick to the UV radiation. 

The rate of momentum deposition from radiation pressure is 
then Prad ~ (1 + TiR)Lincident/c whcrc the term 1 + tir accounts for 
the fact that most of the initial optical/UV radiation is absorbed and 
re-radiated in the IR; tir = Sgas i^m is the optical depth in the IR, 
which allows for the fact that the momentum is boosted by multiple 
scatterings in optically thick regions. Here Sgas is calculated self- 
consistently as the average surface density of the identified clump, 
with AviR ^ 5 (Z/Zq) g~^ cm^ approximately constant over the rel- 
evant physical range of dust temperatures. The imparted accelera- 
tion is directed along the flux vector. In Paper I we discuss numer- 
ous technical aspects of this implementation - such as the effects of 
resolution, photon leakage, and how the momentum is discretized 
- and show that essentially all our conclusions are robust to uncer- 
tainties in these choices. 

The direct momentum of SNe ejecta and stellar winds PsNe 
and Pw are similarly tabulated from STARBURST99 and injected 
as an appropriate function of age and metallicity to the gas within 
a smoothing length of each star. This source of turbulent energy is 
almost always smaller than that due to radiation pressure discussed 
above. In some cases, however, in particular in dwarf galaxies, the 
work done by bubbles of gas ^hock-heated by supernovae and/or 
stellar winds is dynamically important; this is discussed below. 

(2) Supernova and Stellar Wind Shock-Heating: The gas 
shocked by supernovae and stellar winds can be heated to high tem- 
peratures, generating bubbles and filaments of hot gas. We tabulate 
the Type-I and Type-II SNe rates from'M annucci et ar] p006') and 
STARBURST99, respectively, as a function of age and metallicity 
for all star particles, and stochastically determine at each timestep 
if a SNe occurs. That is, the SNe are resolved discretely in time (as 
opposed to continuous energy injection). For each SNe, the appro- 
priate thermal energy is injected into the gas within a smoothing 
length of the star particle. Similarly, stellar winds are assumed to 
shock locally and so we inject the appropriate tabulated mechanical 
power L(t,Z) as a continuous function of age and metallicity into 
the gas within a smoothing length of the star particles. The spe- 
cific energy of these stellar winds is large for young stellar popula- 
tions in which fast winds from massive stars dominate, but declines 
rapidly at later times when slower, AGB winds dominate. 

(3) Gas Recycling: Gas mass is returned continuously to the 
ISM from stellar evolution, at a rate tabulated from SNe and stellar 
mass loss in STARBURST99. The integrated mass fraction recy- 
cled is ^ 0.3. 

(4) Photo-Heating of HH Regions: For each star particle, we 
tabulate the rate of production of ionizing photons from STAR- 
BURST99. Starting from the nearest gas particle and working our 
way outwards radially from the star, we then ionize each gas par- 
ticle which is not already ionized until the photon budget is ex- 
hausted (using the gas and stellar properties to determine the pho- 
ton production rate needed to maintain the particle as fully ionized). 
Gas which is ionized is immediately heated to ~ 10"^ K, unless is 
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already above this temperature; moreover, the gas is not allowed to 
cool below 10^ K until it is no longer in an HII region. This method 
allows for overlapping, non-spherical HII regions that can extend 
to radii ~kpc. 

(5) Long-Range Radiation Pressure: Radiation pressure 
from photons absorbed in the immediate vicinity of stars is cap- 
tured in mechanism (1). However, photons that escape these re- 
gions can still be absorbed at larger radii. For each star particle, we 
construct the intrinsic SED (Ljy) as a function of age and metallic- 
ity; we then use the local density and density gradients to estimate 
the integrated column density and attenuation of the SED using 
p(hsm\ + |Vlnp| ^), where hsmi is the smoothing 
length and is the frequency-dependent opacity (assuming dust 
opacities that scale with metallicity, as in (1) above). The resulting 
"escaped" SED gives a frequency-dependent flux Fu that is prop- 
agated to large distances. We construct a force tree for this long 
range force in an identical fashion to the gravity tree, since after 
attenuating the flux near the star particle, the stellar flux is assumed 
to decrease oc r~^. Each gas particle then sees an incident net flux 
vector F'j^, integrated over all stars in the galaxy. 

Extensive numerical tests of the feedback models are pre- 
sented in Paper I and Paper II. In the Appendix, we present addi- 
tional numerical and convergence studies specific to galactic winds. 



3 WIND MORPHOLOGIES & PHASE STRUCTURE 
3.1 Morphologies with Feedback 

Figure [T] shows the stellar and gas morphology for each galaxy 
model in our standard simulations with all feedback mechanisms 
included. All cases exhibit clear outflows. The stellar maps show 
a mock ugr composite image in which the spectrum of each star 
particle is calculated from its age and metallicity, and the dust ex- 
tinction/reddening is accounted for from the line-of-sight dust mass 
(following [Hopkins et al.|200 5 ). The gas maps show the projected 
gas density (intensity) and temperature (color, with blue represent- 
ing cold molecular gas at T < lOOOK, pink representing the warm 
ionized gas at ^ 10"^ — 10^ K, and yellow representing the hot, X- 
ray emitting gas at > 10^ K). 

The HiZ models result in disks with massive (> 10^ Mq) 
~kpc-scale clumps, similar to "clump-chain" systems. Violent out- 
flows arise from throughout the disk driven by the very high rate 
of star formation, with M* > lOOMQyr"^ The outflow is clearly 
multi-phase: there is a volume-filling hot component (largely from 
hot gas "venting" out of SNe-heated regions and from "fast" OB 
stellar winds), but most of the mass is in dense clumps, filaments, 
and streams. We quantify this more in ^3.5| 

The Sbc models also produce clumpy gas structure but on 
smaller scales and with the global structure more closely trac- 
ing spiral arms. This is similar to observed dwarf starbursts (e.g. 
NGC 1569 or 1313). Again, a multi-phase super- wind is evident, 
arising from the reasonably high specific star formation rate of 
M*/M* -0.2-1 Gyr"^ 

In the MW models, which are much more gas poor than the 
HiZ and Sbc models, the gas in the edge-on image traces the stel- 
lar spiral and bar instabilities, albeit with filamentary and clumpy 
GMC structure internal to the arms. The outflows here are less 
clearly multi-phase; with the low specific SFR and gas fraction, 
much of the material is directly "vented" hot gas, rather than accel- 
erated cold material. 

The SMC models have a patchy ISM typical of dwarf irregular 



systems, with large overlapping SNe bubbles. Outflows above the 
thick disk appear as a mix of venting hot gas and shells of warm 
material accelerated by the hot gas; the latter appear visually as 
loops, shells, and filaments of gas. 

3.2 Effects of Each Feedback Mechanism in Hirn 

Figure |2] shows the edge-on wind morphologies of models that in- 
clude different feedback mechanisms, but are otherwise identical. 

First consider the HiZ model. If we view the disk face-on (see 
Paper II), the clumpy morphology seen in Figure[T]is evident in all 
cases (this arises from violent gravitational instability); however, 
without momentum from radiation pressure, the characteristic sizes 
of individual clumps collapse from ^kpc to ^-^pc and the gas piles 
up at densities ^ lO^cm"^. Supernovae, stellar wind momentum, 
and HII photoionization pressure are not sufficiently strong to resist 
the collapse of clumps beyond a critical threshold and the cooling 
time in these dense regions is much too short for thermal pressure to 
be important. Some hot gas is produced from star clusters that have 
nearly exhausted their gas supply. This vents directly out of the disk 
producing the hot gas seen in the right-most image in Figure [2] but 
the warm/cold material in the wind - most of the wind mass - is 
almost totally absent. 

If we instead turn off sources of gas heating (SNe, stellar wind 
shock-heating, and HII photo-heating), the primary difference rel- 
ative to our standard run is in the volume-filling diffuse phase (the 
properties interior to dense clumps are nearly identical). The diffuse 
gas is much cooler (note the absence of yellow gas in the middle 
panel), but the warm clumps in the outflow - most of the wind mass 
- remain, directly accelerated out of the disk by radiation pressure. 

For the MW-analogue model without radiation pressure, 
GMCs again collapse to significantly smaller sizes, although the 
runaway is not as severe because "fast" stellar winds and HII heat- 
ing are able to somewhat suppress runaway GMC collapse (see Pa- 
per II). Some wind material in the cold phase is also absent, but 
this is not the dominant wind phase in the full simulations, so the 
net difference is much weaker than in the HiZ case. Removing the 
sources of heating has a comparable effect on the SFR as removing 
the radiation pressure (Paper II), but it removes the volume-filling 
"hot" gas and dramatically suppresses the wind (middle panel in 
Fig. [2]); there is thus little material accelerated by radiation pres- 
sure alone in MW-like models. 

Figure|2]shows that the Sbc case lies, as expected, somewhere 
between the MW and HiZ cases in terms of the contribution of 
the different feedback processes to the galaxy-scale outflow. In the 
SMC-like model, on the other hand, the average gas density is very 
low (< 0.1 cm~^), approaching the regime in which the cooling 
time can be comparable to the dynamical time. In addition, with the 
lower densities and metallicities the IR optical depths are not large. 
As a result, with only radiation pressure present, there is no hot gas 
and essentially no wind. Turning off SNe heating, by contrast, has 
a large effect - the volume filling factor of hot gas "bubbles" drops, 
and the wind is basically absent. 

3.3 Appearance in X-Rays 

X-ray observations of canonical systems such as M82 provide a 
strong probe of the shock-heated phase of galactic winds. We there- 
fore consider some of the X-ray properties of our simulated galactic 
outflows. For convenience, rather than making a detailed mock ob- 
servation corresponding to a given instrument, sensitivity, redshift. 



© 0000 RAS, MNRAS 000, 000-000 



Stellar Feedback and Galactic Winds 5 




3 kpc 3 kpc 3 kpc 3 kpc 



Figure 1. Morphology of the gas & stars in our standard simulations with all feedback mechanisms enabled; we show face-on and edge-on projections. The 
time in each is ~ 2 orbital times after the simulation begins, when the disk is in a feedback-regulated steady state. Top: Stars. We show a mock ugr (SDSS-band) 
composite, with stellar spectra calculated from their known ages and metallicities and including dust attenuation. Brightness is logarithmically scaled with a 
^ 2dex stretch. Bottom: Gas. Brightness encodes projected gas density (logarithmically scaled with a ~ 4dex stretch); color encodes gas temperature with 
blue being T < lOOOK molecular gas, pink ~ 10"^ - 10^ K warm ionized gas, and yellow > 10^ K hot gas. Each column corresponds to one of our four galaxy 
models. HiZ: A massive, z ~ 2 — 4 starburst disk with > IOOMq yr~^ ; gravitational collapse forms kpc-scale complexes and a clumpy morphology, with 
violent outflows containing dense, cold gas driven by the massive starburst. Sbc: A z ~ dwarf starburst galaxy; the disk is clumpy but more contained in 
global spiral structure, with outflows from the central few kpc producing a multi-phase clumpy wind. MW: A MW-analogue; the gas morphology more closely 
follows stars in a global spiral pattern; the "wind" here is primarily hot SNe-driven bubbles/holes "venting" out, rather than the clumpy filaments/streamers 
seen in the HiZ/Sbc cases. SMC: An isolated SMC-analogue dwarf; the disk is thick with irregular star formation and large bubbles from overlapping SNe; the 
wind is prominent and contains a mix of hot gas and entrained cold/warm material in filaments/loops/arcs. The "smooth" outer disk is ionized and low-density, 
sufficient to prevent collapse at these masses. 

and energy range, we instead quantify the thermal bremsstrahlung where Tgas is the gas temperature and Ue and m are the number 
emission, for which the emissivity per unit volume is density of electrons and ions. At fixed pressure, ux is typically 

dominated by the lowest temperature gas and so equation [T] need 
not refer specifically to X-ray observations. We quantify the con- 
ux oc T^/^^rierii (1) tributions of gas at different temperatures to the thermal emission 
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standard HiZ No Heating (SNe/Wind/HII) No Radiation Pressure 




Figure 2. Edge-on gas morphologies (as Fig.QJ, with different feedback mechanisms enabled. Left: Standard model (all feedback on). Middle: No heating: 
energy injection from SNe, shocked stellar winds, and HII photo-heating is disabled (gas recycling from these mechanisms remains, but is "cold"). The "hot" 
volume-filling ISM is greatly diminished but in the HiZ/Sbc cases the global morphology and much of the wind mass (in the cold/clumpy phase) remains. In 
the MW/SMC cases, removing heating eliminates most of the wind mass (including warm gas, previously accelerated by hot gas). Right: No radiation pressure 
momentum flux (local or long-range). The hot medium remains, but in the HiZ/Sbc cases most of the cold/warm clumpy material in the wind is absent (hot 
gas simply vents). The wind in the MW/SMC cases is less strongly affected, but within the disk, GMCs collapse much further than they would with radiation 
pressure, giving a higher SFR at the same absolute wind mass. 
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MW model is less clearly wind-like but is instead consistent with a 
diffuse halo. 

The top panel of Figure|4]shows dLx /din Tgas, the contribution 
of gas at different temperatures to the total thermal bremsstrahlung 
emission. This weighting favors dense cool gas given roughly com- 
parable pressures at different temperatures, but there is a rapid cut- 
off below ^ 10"^ K because the gas becomes neutral. Figure|4]shows 
that much of the "x-ray" emission would be contributed by gas at 
^ 10^ — 10^ K with the exact range of temperatures depending on 
the system. Figure |4] also divides dLx /dlnTgas into intervals in ra- 
dial outflow velocity. The material at low Vr < lOkms"^ (includ- 
ing negative Vr) is generally non-outflowing gas at the cooler end 
of the temperature distribution. At higher radial velocities, the ther- 
mal emission tends to come preferentially from higher-temperature 
gas. We show below that in the very high- velocity wind, much of 
the mass is in the form of warm gas, but the thermal emission shown 
here is primarily sensitive to the hotter material. 

We caution that the dominant cooling emission for gas with 
Tgas - 10^- 10^ K probably comes from metal lines, not thermal 
bremsstrahlung. Our treatment of enrichment in the current models 
is relatively simple and just traces the total metallicity (improved 
models will be presented in future work in preparation); therefore it 
is difficult to make detailed predictions for the line emission. How- 
ever, we can use the compiled tables presented in |Wiersma et al.| 
( |2009| ), who use a suite of CLOUDY calculations to determine the 
total metal-cooling luminosity as a function of He, Z (assuming 
solar abundance ratios), Tgas and the photo-ionizing UV background 
(for simplicity the z = background is used), to calculate the metal 
cooling luminosity. We plot the morphology of this emission as 
well in Figure[3] It is broadly similar to the thermal bremsstrahlung 
emission, as expected; however, the metal-line emission is signif- 
icantly more clumpy (and in the Sbc case more concentrated in 
vertical filamentary structures as well) - closer to what is actually 
observed in systems like M82. This is because the metals are in- 
homogeneous. We caution, however, that metal diffusion and other 
small-scale mixing processes are not followed here, so the degree 
of inhomogeneity may be over-estimated. 



Figure 3. Left: Galactic wind thermal emission morphology, as a proxy for 
X-ray emission. The maps show edge-on images of our standard models for 
each disk (with all feedback enabled). Intensity here encodes the projected 
bremsstrahlung emissivity (eq.QJ. A roughly bi-conical, clumpy/multiphase 
wind is typical, similar to that observed in, e.g., M82. Right: Same, but now 
showing the metal cooling luminosity. Note the enhanced dumpiness from 
the inhomogenous metal distribution. 



below. To start, however. Figure [3] shows the observed morphol- 
ogy of the winds for our fiducial simulation of each galaxy model, 
with surface brightness weighted (logarithmically) using the ther- 
mal emission in equation [T] (using the known temperature, density, 
and ionization information in the simulations). Very broadly, the 
appearance of the winds in Figure [3] is similar to that seen in Fig- 
ure[T] but the clumpy structure of the winds is more evident and the 
global morphology is more obviously bi-conical (since much of 
the in-plane material is neutral and so does not appear here). Over- 
all, the morphologies of the Sbc and SMC-like models are quali- 
tatively similar to observations of systems like M82. The HiZ case 
is more spatially extended, making the outflow less columnar. The 



3.4 Velocity and Column Density Distributions 

Figure [4] (middle panel) also plots the radial velocity distribution 
of gas in the simulations, specifically the mass per unit radial ve- 
locity dm/dvr. There is a concentration of gas at small |vr| which 
is the rotationally- supported disk (this includes almost all the star- 
forming gas). This "spike" is symmetric about Vr — 0; the winds 
appear as the large excess of material at Vr ^ 0, extending to 
Vr > 1000 km s~^ even in the SMC-like case. 

In each galaxy model, the winds span a wide range in outflow 
velocity - there is no single "wind velocity" (contrary to what is 
assumed in many models). The distribution of mass in dm/dvr has 
a broad plateau up to some "turnover" velocity, above which it falls 
off more rapidly; we can approximate it with 

drrimnd 1 

dVr ^ 1 + (V,/V0)^ 

with a - 3 - 7 and vo - 800, 500, 200, and 200 kms"^ in the HiZ, 
Sbc, MW, and SMC cases, respectively. For comparison, the max- 
imum circular velocities of these systems are Vmax ~ 230, 86, 190 
and 46 kms~^ In the starburst systems where radiation pressure 
is critical for much of the wind (HiZ/Sbc), the turnover velocity is 
a few times Vmax (comparable to the escape velocity) - this is ex- 
pected for radiatively accelerated winds. In contrast, the maximum 
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Figure 4. Top: Thermal bremsstrahlung emission- weighted distribution of gas temperature dLx /JlogTgas. For each galaxy model (columns), we show results 
for our standard simulation with all feedback mechanisms (time-averaged over the simulation). Different lines in a given panel correspond to gas with different 
radial velocities (relative to the center of mass). The emission from more rapidly outflowing wind material tends to be dominated by the hotter gas with 
T 10^-8 K. Middle: Mass -weighted distribution of gas outflow radial velocities dm/dvy. The distribution for all gas is divided into the contribution from 
different phases: cold (primarily molecular) {T < 2000 /T), warm (primarily ionized) (2000 <T <4xlO^K) and hot diffuse (T > 4x 10^ K) gas. The high 
Vr part of the distribution is the wind material, with a distribution that can be fit by equation [2] This outflowing gas consists primarily of a mix of warm and 
hot gas with some colder (likely molecular) material. The warm material usually dominates the mass, and is in the form of the dense clouds, filaments, and 
shells seen in Figure^ Bottom: Mass- weighted density distribution for the unbound wind material, divided into different phases as in the middle panel. The 
cold/warm/hot gas dominates at high/intermediate/low densities, respectively. The "warm" material at very low densities, <C lO^'^cm"^ is previously "hot" 
material that has adiabatically cooled as it expands; these low densities arise artificially because we do not include an IGM into which the wind expands. 



velocities in the lower-SFR systems appear to asymptote to a max- 
imum at a few hundred kms~^ - such a roughly constant velocity 
is, in turn, the expectation for thermally heated galactic outflows. 

Observationally, the column density distribution of the gas is 
probed directly by absorption line measurements (e.g., Na and Mg 
lines). Figure [5] quantifies this distribution in our simulations. To 
calculate this, we determined the average column density distribu- 
tion calculated over ^100 different projections of the galaxy, each 
time calculating the integrated line-of-sight column along ^10^ 
sightlines through an aperture of fixed physical radius about the 
galaxy center (uniformly sampling the aperture area). Note this 
calculation is simply the total hydrogen column density - we do 
not include the ionization corrections that are critical to comput- 
ing realistic line strengths. The very large columns in Figure [5] 
come from the few percent of sightlines which go through the 
disks, but the bulk of the column density distribution is domi- 



nated by the winds. Figure [5] shows that there is a broad range of 
columns with a roughly power-law behavior in Nh at low values 
(expected for an expanding, constant M wind); the characteristic 
columns range from - 10^^ - 10^° cm around dwarf galaxies to 
^ 10^^ — lO^^cm"^ around massive systems. 

3.5 Cold, Warm, and Hot Outflows 

The mass distribution as a function of radial velocity in Figure |4] 
divides the gas into three traditional phases with a simple temper- 
ature cut: "cold atomic-i-molecular" gas (T < 2000^), "warm ion- 
ized" gas (2000 <T <4xlO^K) and "hot" gas {T >4x 10^ K). 
The wind is a mix of warm gas (containing most of the mass) and 
volume-filling hot gas, with a few percent contribution from cold 
gas at all velocities. This is shown directly in the lower panel of 
Figure |4] where we consider the density distribution of wind mate- 
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Figure 5. Cumulative distribution of gas column density (fraction > Nh) 
along random sightlines within different radii around each galaxy in our 
simulations with all feedback mechanisms. The results are averaged over 
time and projected disk angle, for random sightlines within a projected cir- 
cular aperture of physical radius r around the disk center (labeled). The few 
percent of sightlines at large columns (> lO^^cm"^) go through the disk 
itself; the remainder are dominated by wind material. Wind material with 
Nh ^ 10^^ cm~^ has a large covering factor within a few tens of kpc, and 
there is a broad range of columns ~ 10^^ — 10^^ cm~^. 



rial (defined as all unbound gas, see below), divided into these same 
phases. The wind includes material with a broad range of densities 
cm ^, with high-density material at ;i > 1 cm ^ pri- 
marily "cold," intermediate-density material (n ^ 0.01 — 1 cm~^) 
"warm," and low-density material (n < 0.01 cm~^) "hot." 

The instantaneous phase structure in the wind is not the same, 
however, as the phase structure of the material at the time it was 
initially accelerated into the wind. Its multi-phase origin is evi- 
dent in the fact that the warm gas is highly inhomogenous; most 
of the mass is at a density ^ 20— 100 times higher than would 
be expected if it were smoothly distributed over a spherical shell 
(the formal clumping factor is (p^) / (p)^ ~ 100 in each case). 
Some of the warm gas - especially in the larger and more dif- 
fuse loops and shells, is material from the diffuse "warm-phase" 
of the disk disk (which has both sizeable mass and volume filling 
factors) directly accelerated out of the disk. But the large tail of 
very low-density (n < 10~^cm~^) volume-filling warm gas is ma- 
terial initially launched as "hot" gas, which cools adiabatically as 
it expands. And much of the dense clumpy /filamentary structure is 
material initially launched as "cold" gas, which rapidly becomes 
photo-ionized as it expands and its average density falls. 

There are several important numerical caveats to the micro- 
scale phase structure in the wind shown here. In particular, accel- 
eration of cold gas directly by hot gas is difficult to treat accurately 
in SPH codes, and comparison with adaptive mesh schemes sug- 
gests that we are likely to under-estimate the true amount of cold 
material in the wind (Spring el|2010[|Keres et al.f 2011). Similarly, 
we do not capture the thermal evaporation of cold gas to form hot 
gas or the condensation of hot gas to form cool gas. And the adi- 
abatic expansion of the gas to extremely low densities is possible 
because our simulations do not include a realistic IGM into which 
the gas would expand. As a result, the phase distributions shown in 
Figure |4] are unlikely to be quantitatively accurate. 



4 MASS-LOADING OF GALACTIC WINDS 

4.1 Effects of Each Feedback Mechanism in Turn 

From the broader perspective of galaxy formation, the most im- 
portant integral property of galaxy winds is probably the total wind 
mass outflow rate. This is typically quantified in terms of its ratio to 
the star formation rate Mwind/M*, i.e. the wind mass per unit mass 
of stars formed. Figure |6] shows the wind masses of our models, in 
a series of otherwise identical simulations with different feedback 
mechanisms enabled or disabled in turn. This quantifies the impor- 
tance of a given feedback process in driving the wind in our simula- 
tions. We compare the total mass in the wind Mwind as a function of 
time, relative to the integrated mass in stars formed since the start 
of the simulation Mnew - the ratio of these quantities Mwind /Mnew 
defines the average "global" wind efficiency (Mwind) /(M*). If we 
evaluate the instantaneous efficiency or average in narrow time 
bins, we obtain similar results but with increasing scatter. 

In these calculations, we define the "wind" as material with a 
positive Bernoulli parameter b= (v^ + 3 — Vesc)/2, i.e., material 
that would escape in the absence of additional forces or cooling]^ 

Figure|6] shows that for the HiZ galaxy model, the photon mo- 
mentum plays a critical role in driving much of the wind mass. By 
contrast, removing the sources of hot gas (SNe, stellar winds, and 
HII heating) has little effect on the total wind mass. With all feed- 
back mechanisms operating, a multi-phase wind is blown with an 
outflow rate between ^0.2—1 times the SFR. Without local and 
long-range acceleration by photon momentum, however, the wind 
efficiency drops by a factor of ~ 5, both because the absolute wind 
mass drops and because the star formation rate increases signfi- 
cantly. The long-range acceleration due to the ambient radiation 
field is especially critical in driving the wind: we show in Paper II 
that a moderate increase in the UV escape fraction from the galaxy 
can increase the wind mass by factors of ^ 2 — 5. By contrast, re- 
moving the local IR radiation pressure support within star forming 
regions has a much weaker effect on the wind (though we have 
shown in Paper I that it is critical for the support of GMC com- 
plexes themselves). The winds in the HiZ model therefore seem 
not to be directly launched from GMCs, but rather are the result of 
a two step process (Murra y et al.|2011 ). First, UV and IR radiation 
pressure lofts up gas parcels. Second, these parcels are then contin- 
uously accelerated by photons emerging from the entire disk, rather 
than only those from their natal GMC. 

For the MW model, our fiducial calculations produce a wind 
with an outflow rate a couple times the SFR|^The relative role of 
the hot wind phase is much more prominent than in the HiZ case. If 
we remove all of the thermal feedback (SNe, stellar winds, and HII 



^ We have also calculated the wind mass using a more observationally ac- 
cessible definition: material that has reached at least > 500 pc above the 
disk, with absolute outflow velocity > lOOkms"^ This gas may or may 
not all be formally unbound, and the exact "cuts" in z and v are arbitrary, 
but we find that this definition of wind mass is generally within ~ 20% of 
the wind mass defined with the Bernoulli parameter. Moreover, the results 
are essentially unchanged for moderate (factor ~ 2) changes in the thresh- 
old values of z and v. 

^ The high wind mass-loading at early times (Mwind ~ lOM*) is partially 
an artifact of the initial conditions. This occurs in part because it takes 
~ 2 X 10^ yr for the SFR to reach equilibrium, and partly because the pre- 
simulation bulge-i-disk stars are initialized with a broad age range, so a small 
fraction are sufficiently young to have a large feedback effect before the 
SFR rises much (and because the galaxy is so gas-poor, this has as a pro- 
portionally large effect). 
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Figure 6. Galactic wind mass-loading efficiency (= Mwind/^new; where Mwind = / ^wind and Mnew = / M^) for each galaxy model. We compare our 
standard simulations with all feedback mechanisms to simulations with different feedback mechanisms enabled in turn. Models with no feedback produce 
no measurable wind mass. Since SNe dominate the production of wind material by hot gas (relative to HII heating and stellar winds), the "Rad Mom Only" 
models shown here are very similar to "no SNe" models. Top Left: HiZ model. Net outflow rate is ~ M* . Most of the wind comes from "cold" acceleration by 
radiation pressure, and is present even with no SNe/stellar wind/HII heating terms (though these do contribute a non-negligible fraction from "venting," seen 
in the no radiation pressure case). Top Right: Sbc model, with Mwind ~ 3 — 5M* . Heating mechanisms contribute comparably to radiation pressure in driving 
the wind. Note that the "standard" model wind mass is significantly larger than the sum of the wind mass with just heating and just radiation pressure enabled 
separately - there is a strong non-linear coupling between the different feedback mechanisms that enhances the strength of the wind. Bottom Left: MW-like 
model, with Mwind ~ 1 — lOM*. Here the hot gas produced by SNe dominates the wind driving. Bottom Right: SMC-like model, with M^md ~ 10 — 20 M*. 
Hot gas is again critical, with SNe even more dominant. 



regions), there is essentially no wind! And of the heating mecha- 
nisms, SNe have the largest effect on the wind properties. Absent 
thermal heating, the SFR is still regulated at nearly the same level 
by the radiation pressure feedback (as in Paper I), but material is 
not being driven into a super-wind. Radiation pressure thus plays 
an important role in self-regulating star formation in the disk, but 
is not, unlike in the HiZ models, strong enough to directly accel- 
erate much gas far out of the disk. Rather, radiation pressure lofts 
material out of GMCs, after which the gas is accelerated primarily 
by the hot gas pressure into a galaxy-scale outflow. This contrasts 
with the HiZ case, where the acceleration above the disk results 
from radiation pressure rather than from hot gas ram pressure. 

The Sbc case, as before, lies somewhere in between the HiZ 
and MW cases. The net efficiency is a few times the SFR, with ra- 
diation pressure and hot gas contributing comparably to the outflow 
mass. Note also that Figure[6]clearly shows that the total wind mass 
in the case with all feedback enabled ("standard") is significantly 
larger than the sum of the wind mass with hot gas alone or radiation 
pressure alone. The full wind efficiency stems from the non-linear 
interaction between different feedback mechanisms. 



For our SMC-like dwarf galaxy model, the simulation with 
all feedback mechanisms included drives an outflow of ^ 10 — 20 
times the SFR (Fig.[6|. Just as in the MW model, SNe are the most 
important driver of the galactic outflow. In particular if we remove 
all sources of hot gas (SNe and stellar winds), we again see es- 
sentially no significant super-wind. By contrast, without any feed- 
back from radiation pressure, the wind mass is reduced by a small 
amount (^ 30 — 40%) and the average SFR increases by a similar 
factor, so the net outflow efficiency is only a factor ^1.7 lower. 



4.2 Wind Efficiency versus Galaxy Mass 

Figure [7] summarizes the efficiency of the winds in our standard 
simulations as a function of global galaxy properties. We plot the 
wind efficiency at several different times (evaluated every 10^ yr), 
for each standard disk model, as a function of several parameters: 
the instantaneous SFR, the maximum circular velocity of the disk 
Vmax, and the stellar and baryonic galaxy masses. The order-of- 
magnitude scatter in wind efficiencies even at fixed galaxy proper- 
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Figure 7. Summary of the efficiency of stellar winds, defined as the wind mass loss rate divided by the SFR, Mwind/^* • Each of our galaxy models is shown, 
for simulations with all feedback mechanisms included. The instantaneous efficiencies (the derivative of Fig.[6j are shown at a number of times (black); the 
time-averaged value and ±1 cr scatter over the duration of the simulation is also shown (orange with error bars). The best-fit power-law relation is also plotted in 
each panel (dotted line). Top Left: Outflow rate versus SFR. The best-fit is sub-linear: Mwipd,oc ^O-^^io.os (^^q j^^p j^ight: Efficiency versus the maximum 



circular velocity of the galaxy rotation curve. Our best-fit is Mwind/^* oc V^^x^^'^ (eq.|4|. For comparison, we also show the scaling from cosmological 

simulations that are adjusted to match high-redshift IGM enrichment and approximately match the observed galaxy mass function jOppenheimer & Davej 
u^rx^ii^ 1 • ""TTrTTTr, ifil / ,v /,v \ /T. _i N_i ^ 1 1 , rT.1 • • . -,1 • ■ , — —r — Tr~. T~r^ 



2006 



Oppenheimer et al. 2010)0 (Mwind/^*) = (Vmax/300kms ^) ^red dashed). This is consistent within 1 a with our numerical results. Bottom Left: 



Wind efficiency versus galaxy baryonic mass. Bottom Right: Wind efficiency versus galaxy stellar mass. 



ties is clear, but there are also clear trends in the average efficiencies 
with each galaxy property. 

The wind mass loading is strongly correlated with the galaxy 
SFR, but with a sub-linear power. Fitting a single power-law in 
Figure [TJgives 

^"--3^0^ (m^) 

The sub-linear behavior is particularly interesting, as is implies a 
lower efficiency at higher absolute SFR: Mwind/M* oc M~^^. Since 
observations find a SFR- stellar mass relation M* oc M^^ in star- 
forming systems ( Noeske et al.|2007) , this is equivalent to a lower 
efficiency at higher masses 

This is evident when we directly plot the efficiency versus 
baryonic mass, stellar mass, or maximum circular velocity. The 
strongest correlation appears to be between efficiency and Vmax- 
This has the form 

\ M* / ^ ^VlOOkms-W ^ ^ 

This is in reasonable agreement with the expectation from simple 



dimensional arguments based for a momentum-driven wind with 
a characteristic velocity of Vmax ( [Murray et al.||2005] l. If the mo- 
mentum flux in the outflow scales with the luminosity of massive 
stars (oc L oc M*), then the outward momentum Mwind V'esc oc M*, so 
Mwind/M* oc V~sii- Note that this scaling can also apply if ram pres- 
sure due to hot gas is the primary mechanism driving mass (i.e., 
warm/cold gas) out of the galaxyjjThere is a small, but real devia- 
tion of the HiZ model from the extrapolation of this scaling, which 
is a consequence of an additional residual dependence on Sgas that 
we discuss below. 

The scaling of Mwind/M* with Kiax in Figure [t] and equa- 
tion |4] is in plausible agreement with the scaling suggested in 
cosmological simulations in order to produce a good match be- 
tween the predicted and observed galaxy stellar (and baryonic) 
mass functions as well as high-redshift IGM enrichment patterns. 
jOppenheimer & Dave| p006"] ) and |Oppenheimer et al.|p010| (see 



^ We should note that, although it is not obvious because of the large scatter 
in instantaneous mass-loading rates, the best-fit trend in Figure ^is robust 
to the exclusion of any individual simulation. 
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Figure 8. Best-fit scaling of the wind mass rate as a function of galaxy 
properties in our numerical simulations (eq. [sj: roughly M^inci/^* oc 
V^^ Sgas^^ (where the quantities are evaluated at each wind launching ra- 
dius R in the galaxy). The dotted line shows the identical relation. Black 
points show the instantaneous efficiencies (as Fig. [7}, red points the run- 
ning time-average (as Fig. [6] with quantities evaluated at the half-SFR ra- 
dius at each time). The order-of-magnitude scatter in Fig. [7] is reduced to 

— 1/2 — 

factor ~ 2. The scaling Mwind/^* ^ ^gas reflects the fact that in denser 
systems, SNe cool more efficiently and more of the optical/UV radiation is 
absorbed, reducing the efficiency of wind driving. 



also |Dutton & van den Bosch||2009| ) considered a series of pre- 
scriptions for how galactic winds scale with galaxy properties, 
and ultimately found agreement with observations for a scaling 
Mwind /M* = 3 ( Vmax/ 1 00 km s~ ^ ) ~ \ which lies in thebest-fit range 



of our simulations (as shown by the red line in Fig.[7j[jWe stress 
that we have not tuned or adjusted any parameters in our simula- 
tions to produce either the normalization or slope of this relation. 

Qualitatively similar behavior appears when we consider the 
wind efficiency versus baryonic or stellar mass, albeit with a 
slightly less significant correlation. 



Mwind 



-(0.3-0.65) 



-(0.25-0.5) 



(5) 
(6) 



Given the observed (stellar and baryonic) Tully-Fisher relation 
M* oc Vrnax; these correlations are fully consistent with the Vmax 
correlation; likewise given the M* — M* relation discussed above, 
these are consistent with the wind mass-SFR relation. This consis- 
tency is expected, given that the simulated galaxies lie on both the 
observed Tully-Fisher relations (by design) and SFR-stellar mass 
relation. 



^ [Oppenheimer & Dave|j2006) ; |Oppenheimer et al.H2010) actually do not 
resolve Vmax and instead adopt a proxy for the velocity dispersion a based 
on a galaxy mass estimated from a friends-of-friends object finder. With- 
out a matching cosmological simulation, it is difficult to exactly convert 
between this and the quantities measured here, so there is some ambiguity 
in whether the normalization estimated in the quoted values agrees exactly 
with that implemented. The scaling with V^^x^ however, should be more 
robust. 



4.3 Predicting Wind Efficiencies: Towards a More Accurate 
Effective Model 



In this section we use the results of our simulations to assess if 
there is a more accurate predictor of wind properties than equations 
|3|6| - such a predictor would be extremely useful for application in 
cosmological simulations and/or semi- analytic models. 

First, we consider the simulations alone, including all of our 
galaxy models and case studies where we have varied other pa- 
rameters such as the galaxy structural properties and how the SFR 
is tied to molecular gas. We isolate a series of properties that are 
measurable in such models (each as a function of radius): SFR, 
circular velocity, escape velocity, total/bulge/disk/gas/halo masses 
and surface densities, scale-heights, and the radii themselves (most 
other properties can be constructed from combinations of these). 
We then compare to the corresponding instantaneous wind mass 
loading Mwind(^,0- Although all the parameters are correlated, if 
we marginalize over the other variables we find three parameters 
that dominate the trends in Mwind with radius, time, and galaxy 
model: the SFR M*, circular velocit}^^, and gas surface density 
Ugas. Accounting for these, we do not find > 3a dependence on 
any other variables we consider. A maximum-likelihood fit to these 
three variables gives a correlation of the form 

Mwind OCM^V/S^as (7) 

with a ^ 1.0zb0.15, /3 ^ -1.1 ±0.25, md-f^ -0.5 ±0.15. From 
this, the correlation with SFR is most significant (> 5 a), and at oth- 
erwise fixed properties scales linearly (Mwind oc M*); the sub-linear 
behavior in Fig. [7] owes mostly to the next most-significant corre- 
lation, with Vc{R, t) (^ 4cr), which is consistent with the expected 
oc Vc~^ Third most important (but still > 3 a) is a residual correla- 
tion with surface density which dominates many of the variations in 
Mwind /M* within a given model as a function of radius and time. If 
we include the normalization in the fit, we find the following result 
for the mass loading 

Vc{R) ^|-(l + ^2)/ S,asW \ -(0.5+^73) 



Mw 



M* 



ilOryi { 



VlOOkms- 



/ Sgas (R) 

UOMq pc- 



(8) 

where rii,r]2,r]3 incorporate the uncertainty from the fits (771 ^ 
0.7 — 1.5, 772 ^ ±0.3, 773 ^ ±0.15). A comparison between this 
fit and the simulation results is shown in Figure [8] 

Note that the scaling s with star formation rate and circular ve- 
locity in equation [7] are consistent with the simple momentum con- 
servation argument discussed after equation |4] The surface density 
scaling is, however, new. It implies that the wind mass loading is 
lower for higher gas surface density systems, all other things being 
equal. Physically, our interpretation of this is as follows. First, con- 
sider low mass systems where SNe dominate the wind driving. As 
shock-heated SN bubbles expand, the momentum of the swept-up 
material increases in time due to the P-dV work done by the hot 
gas. This enhances the wind driving so long as cooling does not 
sap the energy of the SN bubbles. For higher surface density galax- 
ies, however, radiative cooling of SN remnants is more important, 
suppressing the SNe contribution to the wind. This is qualitatively 
consistent with the scaling in equation [8] 

Consider now the opposite limit of massive, dense systems in 
which radiation pressure is the most important wind driving mech- 
anism in our models, in particular the long-range continuous accel- 
eration due to the diffuse radiation field. The gas at large heights 

^ Because we are considering the winds at different galactic radii 
^wind(^,0' we use Vc(R) instead of Vmax = MAX(Vc[R]). 
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above the midplane is generally not optically thick in the IR, but 
still has T > 1 in the UV (and possibly in the optical, if the wind 
is sufficiently dense). Thus the incident radiation pressure force on 
this gas depends on the escape fraction in the UV/optical /esc,uv-opt 
times L/c, so that Mwind oc /esc,uv-optM* V^^). In the Appendix, we 
explicitly confirm that the wind mass-loading in the HiZ and Sbc 
models indeed scales close to linearly with /esc,uv-opt. The escape 
fraction in turn declines for larger disk gas surface densities, be- 
cause the disk is itself what is absorbing the UV flux. This implies 
a wind mass-loading that decreases for higher gas surface densities, 
again qualitatively consistent with the scaling in equation [8] 



5 DISCUSSION 

We have studied the origin of galaxy- scale outflows using a new 
set of numerical methods for modeling stellar feedback in hydro- 
dynamic galaxy simulations; the feedback processes include radia- 
tion pressure on the scales of star clusters; shock-heating, momen- 
tum injection, and gas recycling via supernovae and stellar winds; 
photo-ionization in (overlapping) HII regions; and radiation pres- 
sure produced by the diffuse interstellar radiation field. Our cal- 
culations use the results of stellar population modeling to self- 
consistently include the time-dependence of these feedback pro- 
cesses. We have explored this physics in the context of isolated 
(non-cosmological) galaxy models that range from those motivated 
by massive z ^ 2 galaxies forming stars at ^ 100 — 300 Mq yr~^ 
to models of SMC-like dwarf galaxies. 

By incorporating multiple feedback processes, we find that nu- 
merical simulations produce galaxy-scale outflows with mass load- 
ing factors up to Mwind ^ lOM*. Moreover, the wind mass-loading 
is the largest in the lowest mass galaxies, scaling approximately as 
~ (although with an additional gas surface density dependence 
that we discuss below). This is, to our knowledge, the first time 
that hydrodynamic simulations have demonstrated that stellar feed- 
back processes can self-consistently generate galactic winds with 
sufficiently large mass-loading and the appropriate galaxy mass- 
dependence to match what has been invoked to explain the shape 
of the the galaxy stellar mass function. 

A general feature of our results is that the interaction between 
multiple feedback processes, acting on different spatial and tempo- 
ral scales, is important for driving realistic galactic winds. In par- 
ticular, some feedback processes prevent runaway collapse within 
GMCs, thus maintaining a diffuse phase of the ISM that is more 
susceptible to feedback. In general a different set of processes then 
accelerates the wind material out of the star-forming disk. 

By turning off individual processes in our numerical simula- 
tions, we can identify the mechanisms responsible for most of the 
wind mass-loading in a given galaxy model. In dense gas, typi- 
cal of high-redshift rapidly star-forming galaxies (our HiZ model) 
and local starburst galaxies (local LIRGs and dense nuclei such as 
Arp 220), cooling times even in the "diffuse" ISM are extremely 
short, and thus energy input into the ISM by SNe and stellar winds 
is rapidly radiated away. This heating can neither halt the run- 
away collapse of massive GMC complexes or accelerate significant 
amounts of material into a wind (since most of the gas is in the 
dense medium). Instead, radiation pressure from trapped IR pho- 
tons in star-forming regions serves to dissociate those regions, and 
"lofts up" gas above the disk - there, the gas "feels" the coherent 
momentum flux from the optical/UV/IR photons escaping the other 
star clusters in the disk, and can be continuously accelerated out of 
the galaxy. Most of the wind mass in our models of these dense 



systems is in cold, dense blobs radiatively accelerated out of the 
galaxy; the small amount of hot gas is volume filling but has little 
mass, and is too tenuous to significantly accelerate the dense mate- 
rial. Under these conditions, we find that the net wind efficiencies 
are typically moderate, Mwind 0.1 — 2M*. 

In lower-gas density systems such as the MW and dwarf star- 
bursts such as M82 (analogous to our Sbc model) {{n) ^ 1 cm~^), 
heating mechanisms and direct momentum injection have compa- 
rable effects. Injection of momentum by radiation pressure, 0-star 
winds, and photo-heating when stars clusters have ages < 10^ yr be- 
gins to dissociate their host GMCs and punch holes that allow SNe 
remnants to expand and escape into the diffuse/hot medium. There, 
cooling times are low, and bubbles can grow, overlap, and do PdV 
work on the entrained warm and cold gas, contributing significantly 
to the wind driving. With all feedback mechanisms included, dwarf 
starbursts generate galactic winds with mass-outflow rates ^2 — 5 
times the SFR. 

In our dwarf galaxy model (motivated by the SMC), and in the 
outer extended disks of more massive systems, the gas densities are 
very low ({n) < 0.1 cm~^). The cooling time can be comparable 
to or longer than the dynamical time so that heating mechanisms 
become very efficient; such systems also have low stellar surface 
densities and gas optical depths, so radiation pressure is less im- 
portant. A combination of momentum from UV photons, stellar 
winds, and warm gas pressure from photo-ionized regions stirs up 
star-forming regions, allowing SNe explosions to expand rapidly 
and generate large, overlapping hot bubbles. This generates a char- 
acteristic dwarf irregular morphology with patchy star-forming re- 
gions and hot SNe bubbles (see also Governato et al.|20T0| . These 
hot bubbles push cold shells and loops into a super-wind as they 
overlap, generating a multi-phase (mostly hot/warm) wind with a 
large mass loss rate Mwind 10 — 20M* . Without SNe, there is still 
some hot gas and a "bubble-like" morphology driven by high ve- 
locity stellar (0-star) winds, but SFRs are a factor > 5 larger and 
(absolute) mass outflow rates a factor of ^ 20 lower. 

We have used our simulations to determine improved prescrip- 
tions for wind mass-lass rates (eq. jSj and velocity distributions 
(eq[2]) for use in cosmological simulations and semi-analytic mod- 
els of galaxy evolution (for which the explicit models here cannot 
yet be implemented): Mwind /M* oc Vc^^Sg^s^^. This fit is similar 
to what has been adopted phenomenologically in the past, but in- 
corporates an additional correction such that the mass-loading de- 
creases with increasing gas surface density oc Sg^s^ (at fixed cir- 
cular velocity). In the radiation-pressure dominated regime, this re- 
flects the fact that a dense disk captures more of the primary stellar 
luminosity very close to stars, downgrading photons into the IR. 
These degraded photons can help drive turbulence in very dense 
regions (with tir^ I) but generally do not contribute to acceler- 
ating material above the disk. In the SNe-dominated regime, the 
surface density scaling reflects the fact that more of the SN energy 
is radiated away in denser systems. 

In addition to introducing variations between galaxies with 
the same Vc, the Sgas dependence will introduce some additional 
systematic scaling with mass and circular velocity. At fixed red- 
shift, total disk surface density tends to increase with mass and Vc, 
but at intermediate masses (M* > IO^^Mq) gas fractions also de- 
cline with mass such that Egas (measured within the star-forming 
disk, which is the important scale for winds) is nearly indepen- 
dent of galaxy mass (e.g. |McGaugh||2005l |Hall et al.|[20TT] ). At 
lower masses, however, where Mgas ^ M*, Egas scales with mass as 
^ Vc~^, steepening the mass-loading by an additional power of Vc. 
This may resolve some of the problems noted in |Dave et al.lpOlT) , 
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where a simpler single power-law scalings tends to over-produce 
the number of low-mass galaxies. 

The scaling we have inferred for the mass loss rate is appli- 
cable to a range of galaxies but may break down in two different 
extremes, which we will explore in future work. In the lowest mass 
dwarf galaxies (with M* < \& Mq and Z < 0.01 Zq), the ioniz- 
ing background alone may be sufficient to suppress most structure, 
metallicities may be sufficiently low as to prevent molecule forma- 
tion except under extreme conditions, and a single SNe can clear 
out the gas supply; our simulations do not apply to these circum- 
stances. At the opposite extreme, in the most dense regions of nu- 
clear starbursts and material surrounding AGN, even the "diffuse" 
medium can be optically thick in the IR and many of our assump- 
tions about the efficiency of radiation pressure and the (in)ability 
of IR photons to drive winds are likely to be modified. 

We also caution that care is needed when mapping the scaling 
here - which assumes the disk structure is marginally resolved and 
so quantities like Vc and H can be evaluated at the effective radius 
of star formation - to models with poor resolution. For example, 
in the'Oppen heimer & Dave| ( |2006 ) cosmological simulations (and 
many analytic models), /?disk and Vc(^disk) are unresolved, so they 
are implicitly assumed to scale in a simple manner with total galaxy 
baryonic mass (and redshift). If we use the observed Tully-Fisher 
relation and gas fraction- stellar mass relations as noted above to 
convert our proposed scaling to a strict power of mass and redshift, 
the scalings are similar at intermediate masses. But it is in general 
non-trivial to link the quantities relevant to the outflows in our sim- 
ulations (e.g., Fig.jSj to global parameters like halo/baryonic mass 
or halo circular velocity. 

The calculations presented in this paper can clearly be im- 
proved in a number of ways in future work. Most notably our 
treatment of the radiation pressure produced by UV, optical, and 
IR photons is quite approximate. A full calculation including the 
scattering and absorption of both the UV and IR photons is tech- 
nically formidable. We have calibrated our treatment of the diffuse 
radiation field in galaxies using such radiative transfer calculations 
(Paper II) but additional improvements to our methods would be 
valuable. Our calculations also do not include cosmic-rays gen- 
erated in SN shocks, which empirically contribute significantly to 
the pressure-support of MW-like spiral galaxies (Bou lares & Cox| 
|1990j ) and may also contribute significantly to driving outflows 
from such galaxies. And we reiterate that even at our highest resolu- 
tion, there are important numerical caveats (and unresolved mixing 
processes) that may apply to the phase structure of the winds. 

Finally, our calculations have focused on idealized, isolated 
galaxy models and do not include a realistic hot coronal gas com- 
ponent, galactic inflow, or an external intergalactic medium. The 
interaction between galactic outflows, hot halo gas, and inflow- 
ing cold streams is an extremely interesting problem that should 
be studied in detail in future work. Because we have not captured 
this interaction, the outflow rates we calculate are best interpreted 
as outflow rates from the star-forming disk, rather than from the 
virial radius of the host dark matter halo. 
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Figure Al. Effects of allowing for photon "leakage" - i.e. assuming that 
the ISM is clumpy and porous on sub-resolution scales. We plot the wind 
mass-loading versus time (as Figure [6} for an otherwise identical standard 
HiZ simulation. Our standard model ("without leakage") takes the local 
column density estimator as-is, and estimates the escape fraction from the 
spectrum of each star with /esc = exp (— r) for each waveband. The "with 
leakage" model assumes that each star is actually surrounded by an (un- 
resolved) distribution of columns following a power-law distribution (see 
Paper I for details) and uses this to calculate /esc - The allowance for leakage 
only makes a difference when the average opacity is large, at which point 
the contribution (in either case) to the global total luminosity in the relevant 
band is small. As a result it has a negligible effect on our results. 
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APPENDIX A: NUMERICAL TESTS 

Extensive numerical tests of the models used here are presented in 
Paper I and Paper 11. We refer to those papers for details, but briefly 
summarize the results of those numerical tests for galactic winds. 
In Paper I, we note that the ISM is probably inhomogeneous 



below our resolution scale, which can allow photons to "leak out" 
from the region around their parent star at a rate higher than the 
nominal exp(— ro) which we use to attenuate the spectrum when 
calculating the radiation pressure effects. This will weaken the 
short-range radiation pressure, but strengthen the long-range ra- 
diation pressure forces (since more short-wavelength photons es- 
cape, to which the "diffuse" medium has a higher optical depth). 
In Paper I and Paper II, we show that the "escape fraction" (and 
hence modification to these forces) is straightforward to calculate 
analytically, if we assume some a priori functional form for the 
sub-grid local density distribution. In Figure [AT[ we show the con- 
sequences for the winds if we replace our usual estimator of local 
attenuation, /esc = exp (—to), with the appropriately modified /esc 
assuming that a gas particle with mean opacity to really represents 
a complete distribution of opacities with that median but a power- 
law functional form and dispersion of ^ 0.6 dex (chosen to match 
the median from higher-resolution simulations of sub-regions). We 
test this for the HiZ model, since this is the model most sensitive 
to the radiation pressure effects. Adopting a lognormal distribution 
gives a nearly identical result, for all dispersions ^0 — 2 dex (see 
Paper I). In Paper I we showed our treatment of the subgrid escape 
fraction had little effect on the SFR. Here, we see that it also has no 
significant effect on the resulting galactic wind. The reason is that 
for low/intermediate to, the escape fractions are large (and hence 
similar) regardless of sub- structure; for high to the escape fractions 
can be very different, but they are both still small and contribute 
negligibly to the total coupled momentum. 

In Paper II, we note that we can isolate the effects of the 
emergent galaxy SED shape on the winds by comparing our full 
model (in which the emergent flux from each star particle in the 
UV/optical/IR is estimated from the stellar age and local attenu- 
ation) to one in which we simply set the SED shape by hand and 
force it to be the same for all particles. Specifically, we calculate the 
bolometric luminosity and incident force on all gas particles in our 
standard fashion, but simply force the values Luv/^boi, Lopt/^boi, 
Lm/Lhoi to equal some chosen values. 

Figure [A2] compares our standard model, for each galaxy type, 
to two choices for this "enforced SED" model. First, we consider 
a model in which the SED shape is chosen for each galaxy model 
to match observations of galaxies with similar stellar mass, SFR, 
and redshift. The observations are discussed in Paper II, but for 
typical observational counterparts to each of our galaxy models, 
(/uv, /opt, /ir) is approximately equal to (0.05, 0.2, 0.75) for the 
HiZ, (0.13, 0.40, 0.47) for the MW, (0.3, 0.3, 0.4) for the SMC, 
and (0.07,0.23,0.7) for the Sbc models. There is substantial scat- 
ter in otherwise similar galaxies and uncertainty in these quantities, 
but they provide an approximate guide. Reassuringly, the result- 
ing wind masses are quite similar to our "standard model," so it 
is likely that the simulations are capturing the key physics of ra- 
diative acceleration despite lacking a complete on-the-fly radiation 
hydrodynamics. This is consistent with the fact that in Paper II we 
showed that the predicted SEDs from our radiation pressure model 
are similar to those observed. 

We also show in Figure |A2| the wind mass-loading if we 
were to assume a uniform SED shape with a a significantly higher 
UV/optical escape fraction. For all models, this serves to substan- 
tially enhance the coupled momentum for material above the disk 
and so enhances the wind mass (and suppresses the SFR, shown 
in Paper II). This confirms two important points discussed in the 
text. First, on large scales (as opposed to e.g. inside GMCs) the gas 
is optically thin in the IR, and so it is the UV/optical photons that 
dominate the long-range radiative acceleration of winds. Second, 
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Figure A2. Comparison of our standard long-range radiation pressure model (where the attenuation and SED shape of each star particle is calculated) to 
one in which we simply force the SED shape to match a specific empirical template. We show the resulting wind mass-loading versus time (as Figure |6j 
for each of the HiZ/Sbc/MW/SMC models. In each, we compare our standard model (black) in which the SED shapes are self-consistently calculated 
on-the-fly, to a model where the SED shape (Lboi is still determined self-consistently) is forced to a constant fixed value chosen to match the observed 
mean SED shape for similar observed galaxies to each model (blue). These choices amount to a relative proportion of Lboi in the UV/optical/IR bands of 
(0.05, 0.15, 0.8), (0.07, 0.23, 0.7), (0.1, 0.4, 0.5), (0.2, 0.3, 0.5) for the HiZ/Sbc/MW/SMC cases. In each case, the self-consistent and empirically fixed 
models give SFRs in reasonable agreement; this is because the emergent SEDs calculated with the full models tend to agree well with the typical observed 
values (see Paper II). However factor < 2 differences do arise, largely because the empirical model does not allow spatial and/or time-dependent variations 
in the emergent SED (for example, very young stars in the galaxy nucleus tend to be much more obscured than older stars in the galaxy outskirts). We also 
compare a model with fixed UV/optical/IR fractions, but much higher UV/optical fractions than actually observed (red): (0.2, 0.3, 0.5) for the HiZ/Sbc/MW 
cases and (0.4, 0.4, 0.2) for the SMC. This artificially boosts the long-range radiation pressure, which leads to a nearly linear corresponding boost in the wind 
efficiencies. 




Figure A3. Resolution tests. We plot the wind mass-loading versus time 
(as Figure |6j. We consider a "standard" (all feedback enabled) HiZ model 
at a series in mass and force resolution (A/gas is the number of gas parti- 
cles in the star-forming disk, which determines the mass resolution; Cgas is 
the force softening). The wind mass tends to increase at higher resolution 
as venting, mixing, and acceleration of filaments/clumps/streams through a 
background medium can be more accurately resolved. But once we reach 
~ 10^ disk particles and ~ lOpc resolution, the winds as well appear to 
converge (changing by < 10% going to higher resolution). Recall, the Jeans 
length in these galaxies is several hundred pc (Jeans mass > 10^ Mq), so 
all of these models formally resolve the Jeans scales. At lower resolution, 
our prescriptions do not have clear physical meaning. 



force softening (spatial resolution)]^ In Paper II the SFRs for this 
study are shown - they converge relatively quickly. In contrast, the 
winds show a clear, albeit not dramatic, resolution dependence. The 
fact that the higher-resolution runs converge to equilibrium more 
rapidly is not important - it is simply a consequence of the fact 
that most rapidly collapsing spatial scale resolved is on a smaller 
timescale as we go to higher resolution. But going from low to in- 
termediate resolution, the equilibrium wind mass loading efficiency 
approximately doubles. At still lower resolution than we show here, 
the wind mass begins to drop exponentially. This is because the 
simulation can no longer resolve the phases of the ISM and "holes" 
for gas to accelerate out along. In addition, at low resolution the 
columns around each star particle (hence extinction of UV/optical 
light) will be over-estimated, since the dumpiness of the ISM is not 
correctly treated. Moreover at low mass resolution it requires very 
concentrated momentum to accelerate even single particles out of 
the disk - this is why the results appear to be more dependent on the 
mass resolution (A/gas), at least in the range we show here, than on 
the force resolution at fixed mass resolution. Finally, at low resolu- 
tion in particular, it is well-known that SPH codes with energy and 
entropy conserving integration schemes produce artificial "surface 
tension" forces that suppress the acceleration gas streams through 
a porous or unstable medium. Once we reach intermediate reso- 
lution, however, the convergence even in the winds becomes quite 
good. From intermediate to standard resolution the wind mass load- 
ing changes by only ^ 10 — 20%, and from standard to ultra-high 
resolution it changes by < 5%. 



the "boost" in the wind mass-loading is approximately linear in the 
UV/optical escape fraction, as predicted in § |4.3| This drives the in- 
verse dependence of the wind mass loading on gas surface density 
(at otherwise fixed properties). 

Al Resolution Tests 

In Figure |A3| we consider a basic spatial and mass resolution 
study of the winds. We consider an otherwise identical "stan- 
dard" HiZ run, with different particle number (mass resolution) and 



^ For technical reasons to ensure quantities such as random number gener- 
ation were identical in these runs, this set was run on the same set of pro- 
cessors and so does not overlap with our standard HiZ run shown earlier; 
but the results from that run are completely consistent with the standard- 
resolution case here. 
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